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The existence of anharmonic localization of lattice vibra- 
tions in a perfect 3-D diatomic ionic crystal is established for 
the rigid-ion model by molecular dynamics simulations. For 
a realistic set of Nal potential parameters, an intrinsic local- 
ized gap mode vibrating in the [f f f] direction is observed for 
fee and zinc blende lattices. An axial elastic distortion is an 
integral feature of this mode which forms more readily for the 
zinc blende than for the fee structure. Molecular dynamics 
simulations verify that in each structure this localized mode 
may be stable for at least 200 cycles. 

It has been proposed and numerically demon- 

strated that a large amplitude vibration in a perfect 
1-D lattice can localize because of anharmonicity. More 
detailed analytical and numerical investigations of classi- 
cal 1-D anharmonic chains made possible by simple eigen- 
value generating recursion relations have revealed a va- 
riety of stable intrinsic localized modes with frequencies 
outside of the plane wave bands ^-|^ . Recently quantum 
mechanical aspects of ILMs (Intrinsic Localized Modes) 
have been considered ||Jl^ . Some progress also has been 
reported for higher dimensional classical crystal lattices 
with simple nearest- neighbor model interactions Jll| , p^ . 
Particularly important has been the recognition that di- 
atomic crystal potentials like the Born-Mayer-Coulomb 
in 1-D produces an intrinsic gap mode (IGM) between 
the optic and acoustic branches instead of an ILM above 
the plane wave spectrum |l^ . The possibility of IGMs 
in 3-D anharmonic lattices with realistic potentials has 
remained elusive. One approach has been to focus on 
crystal surfaces and edges where harmonic localization 
already plays an important role [p^|-p^ . 

In this Letter we demonstrate with molecular dynam- 
ics simulations that, for sufhciently large vibrational am- 
plitude, anharmonicity can stabilize an IGM in a 3-D 
uniform diatomic crystal with rigid ion Nal potential ar- 
ranged in either the fee or zinc blende structure. By 
developing a self-consistent numerical technique for find- 
ing an intrinsic localized mode eigenvector, we have been 
able to show that for a given gap mode amplitude with 
the same potential that the localization is much stronger 
for a Td symmetry site when compared to an Oh one. 

To construct the stationary localized mode eigenvector 
for the nonlinear 3-D diatomic lattice with long range in- 
teractions, we build on techniques which have been used 
successfully to identify 1-D anharmonic modes ||l3|, one 
of which is the rotating wave approximation. A Fourier 
extension of that idea for vibration in a stationary peri- 
odic mode with fundamental frequency, lo, which includes 
the static and second harmonic for the j*'' particle dis- 



placement, Vi{t), is 
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where r-°\ is the distorted equilibrium position of the 

i*'* particle, and the r^"^'s are the time-independent am- 
plitudes of the different harmonics. The force acting 
on the i*'' particle can also be represented by a simi- 
lar Fourier series. Substituting the coordinate and force 
Fourier series into the classical equations of motion and 
equating the terms with different harmonics gives a time- 
independent system of nonlinear equations 
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The Fourier coefficients, F^"^ , are determined in the usual 
way 

The procedure for finding the anharmonic localized 
mode eigenvector relies on first generating a fictitious dy- 
namics for the Fourier amplitudes and then applying a 
version of dynamical simulated annealing [ pT| . The coef- 
ficients Y^^'^ , rf'^ , and rf'^ are now taken to vary with time 
and obey the following fictitious equations of motion: 



' = F - ^ -f n UJ rriiVi ' 
n = 0,l,2; 
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so that dynamical simulated annealing can be used to 
find the equilibrium values for these quantities rl"-*, r^^-*, 
and vf \ when all f["^ (''^"^) ~ ^' ^™ iterative proce- 
dure to obtain the eigenvector is as follows: 

1. Guess an initial shape for the gap mode eigenvec- 
tor. In our case we use the mass-defect eigenvector 
associated with a harmonic gap mode of frequency 
UJ. The rescaled mass-defect gap mode eigenvector 



gives initial amplitudes r^^'' where r. 



a while 



if\ rp-*, and the initial velocities are set to zero. 

2. Make an MD time step by solving Eqs. (3) and up- 
date the quantities, vf^ and the ■* for fixed cen- 
tral particle amplitude a. The classical equations 
of motion are integrated using the "leap-frog" al- 
gorithm llj] with a time step of 1.35 fs. 

3. Apply a form of dynamical simulated anneal- 
ing jl^. During a MD simulation run, the oscil- 
latory values of the kinetic energies of each of these 
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three objects rj"' per particle for the N particles is 
monitored. (Each object vibrates around its equi- 
librium position r^"''). When the total kinetic en- 
ergy of the system of objects passes through its 
maximal value in the MD fictitious time-evolution, 
•(") 

all object's velocities are set to zero. In this 
way we incrementally move the system closer to its 
equilibrium configuration. 

4. After 40 MD steps as described in ^ and check- 
ing for condition ^ the intrinsic gap mode fre- 
quency Lo is updated by solving the single eq. 
Fq^^ + Lu'^mor^^ = for cj with r^^^ — a. 

5. Verify the correctness and stability of the resultant 
IGM eigenvector with a regular MD simulation. 
Repeat |^ through ^ until the lifetime of the mode 
remains unchanged when the resultant eigenvector 
is used as an initial condition for an MD simulation. 
As long as the frequency of the mode remains in the 
gap we find the procedure described here converges 
to a localized eigenvector. No change in the results 
is observed when this procedure is repeated with a 
1/2 time step. 

To investigate the possibility of a 3-D intrinsic gap 
mode in ionic crystals we use the tabulated rigid ion po- 
tential for Nal from Table's 2 and 4 of Ref. ||l9|. The 
potential has the following form 

i^ij) = ^l^ll- + exp [-k {nj - ro^j)] (4) 

^6 ^8 

ij ij 

with Tij the distance between the ions i and j, Zi and Zj 
are ±1 charges, eo is the permittivity of vacuum, roy is 
the sum of the ionic radii, b and k are parameters deter- 
mined by fitting the thermal expansion and isothermal 
compressibility, Cy and dij are respectively the coeffi- 
cients for the dipole-dipole and dipole-quadrupole in- 
teractions. The lattice constant for the fee structure is 
chosen to be a = 6.35A, within 1% of the experimen- 
tal value [20] . The calculated TO frequency is ujto — 
2.51 X 10^"^ rad/s, within 5% of the experimental value 
of 2.39 X 10^^ rad/s The gap between the optic and 
acoustic branches extends from cj+ = 2.41 x 10^'^ rad/s 
to = 1.46 X 10^3 rad/s. 

In carrying out the MD calculations on a 216-ion cube 
with periodic boundary conditions the method of Sang- 
ster and Dixon has been used to evaluate the Ewald 
sum. For this method the cut-off distance in real space 
is approximately half the length, L/2, of the cube (the 
interaction of sixth neighboring shell is counted), and the 
reciprocal lattice is summed with the convergence param- 
eter of 5.6/L. 

The resulting intrinsic gap mode eigenvector for the 
fee Nal crystal is shown in Fig. 1. Panel (a) identifies the 



first harmonic part of the vibrational eigenvector, , 
localized on a central light Na+ ion and its neighboring 
shells with the central Na"*" ion vibrating in the [111] 
direction. The vibrational amplitude of the central Na"*" 
ion to the nn distance, d, r^^^ jd = a/d = 0.244. The gap 
mode frequency Hes 5.0% below the bottom of the optic 
band. The elastic distortion is characterized by the dc 
part of the mode's eigenvector, r^"'', which is shown in 
panel (b) of Fig. 1 

The above procedure has also been performed on a 
larger array of particles to insure that the periodic bound- 
ary conditions do not influence the results. In a 1000 
ion crystal the IGM with amplitude a/d — 0.244 has the 
same eigenvector as shown in Fig. 1 and its frequency dif- 
fers by 1% from that for the 216 ions Nal crystal. This 
frequency difference is associated with a slightly different 
crystal distortion between the two cases. The MD sim- 
ulation test shows that the IGM remains stable in the 
large crystal and its life-time increases slightly. 

In order to investigate the role of point group symme- 
try on the intrinsic gap mode parameters the zinc blende 
structure has also been tested with the same potential. 
(A local minimum of the lattice energy occurs at the 
slightly larger lattice constant, a = 7.00A.) The IGM 
eigenvector for 216 particles is shown in Fig. 2. Again 
the central ion vibrates in the [111] direction. This mode 
has an amplitude to nn distance, /d — a/d ~ 0.116. 
Again the mode vibrational eigenvector, r^^-*, is localized 
on a central light Na"*" ion and its neighboring shells. 
Note that although the relative mode amplitude is only 
one half that of the IGM shown in Fig. 1, its relative fre- 
quency occurs at the same value in the gap (5.0%) due to 
the larger elastic distortion in the zinc blende lattice near 
the mode center. The maximum vibrational amplitude 
and dc distortion for each of the shells of particles shown 
in Figs. 1 and 2 are given in Table |. 



TABLE I. Maximum amplitude in a shell versus shell index 
for Figures 1 and 2. Both the vibrational amplitude and the 
de distortion are given as fractions of the nn distance. 





fee 




zinc blende 




shell 


vib. ampl. 


dc dist. 


vib. ampl. 


dc dist. 





0.2438 





0.1156 


0.0291 


1 


0.0119 


0.0244 


0.0147 


0.0220 


2 


0.0734 


0.0074 


0.0242 


0.0064 


3 


0.0025 


0.0132 


0.0005 


0.0039 


4 


0.0149 


0.0078 


0.0039 


0.0038 



Typical power spectra of the central particle vibration 
is presented in Fig. 3 for the eigenvectors shown in Figs. 1 
and 2 after about 200 vibrations. These spectra reflect 
the stable vibration of the IGM at frequencies close to 
the values predicted by Eq. (||) both for the fee and zinc 
blende structures. Note that the vibrational mode for the 
Oh symmetry site shows a weak third harmonic while the 
IGM for the Td symmetry site shows all harmonics. 
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FIG. 1. The intrinsic gap mode eigenvector in the Nal crys- 
tal with fee structure. The mode hcis a relative amplitude 
a/d = 0.244 and a relative frequency = 0.950. Panel 

(a) shows the first harmonic part of the eigenvector. The ar- 
rows represent 10 x actual amplitudes for clarity. Panel (b) 
shows the dc part of the eigenvector. The arrows representing 
the displacements are expanded 40 x for clarity. Small cir- 
cles, Na"*" ions; large circles, I~ ions. Although not clear from 
the figure the magnitudes of the dc distortion for the nearest 
I" ion shell arc the same, consistent with an axial distortion 
around the [111] axis. 




FIG. 2. The intrinsic gap mode eigenvector for a hypothet- 
ical Nal crystal with zinc blende structure. The mode has 
a relative amplitude a./d = 0.116 and a relative frequency 
u!/uj+ = 0.950. Panel (a) shows the first harmonic part of the 
eigenvector. The arrows represent 20 x actual amplitudes. 
Panel (b) shows the dc part of the eigenvector. The arrows 
representing the displacements are expanded 40 x . Small cir- 
cles, Na"*" ions; large circles, I~ ions. 
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FIG. 3. The power spectrum of the central particle vibra- 
tion for an IGM in an Nal crystal with either fee or zinc blende 
structure. The eigenvectors shown in Figs. 1 and 2 are used as 
the initial conditions for the MD simulations. The solid curve 
(fee lattice) is shifted down by 10 decades from the dashed 
curve (zinc blende) for clarity. The noise is associated with 
truncational errors. 
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cx/d 

FIG. 4. The frequency of an intrinsic gap mode as a func- 
tion of normalized amplitude, ajd. The left set of data are 
for the zinc blende and the right set for the fee lattice. The 
numerical solutions of Eq. (j^ are represented by the solid 
lines for the [111] and by dashed lines for the [110] crystal di- 
rection. The result of MD simulations are given by the open 
diamonds for the [111] direction and by open circles for the 
[110] direction. 



Figure 4 shows the IGM frequency verses the ampli- 
tude for the two structures under investigation. The left 
sets of data are for zinc blende and the right sets are for 
the fee lattice. These results indicate that the sym- 
metry site appears to support more anharmonicity in the 
sense that for a given vibrational amplitude the frequency 
of the IGM drops farther into the forbidden gap and has 
a larger elastic lattice distortion around the IGM center 
[compare Figs. 1(b) and 2(b)]. 

The general behavior of the IGM frequency versus am- 
plitude shown in Fig. 4 is similar for both structures — 
with increasing mode amplitude the frequency drops far- 
ther into the gap. An amplitude threshold is evident. 
An IGM which has its frequency about 5% below the 
bottom of the optics band (corresponding to the middle 
regions of the curves in Fig. 4) has the longest lifetime of 
~200-250 vibrational periods. If its frequency drops far- 
ther into the forbidden gap (~10%) the mode's lifetime 
decreases to '--^100 periods. At the opposite small am- 
plitude limit its hfe-time again decreases (~40 periods) 
presumably because of the strong coupling between the 
IGM and the nearby plane waves. 

In both systems the elastic distortion associated with 
the IGM has lower symmetry than the corresponding 
point group symmetry of the crystal. To recover the full 
crystal point group symmetry for this perfect crystal it 
must be possible with a different set of initial conditions 
to rotate the IGM about the equilibrium lattice site. To 
test this idea we have examined IGM excitations along 
the three different crystal directions. Although a stable 
IGM mode does not appear for initial excitation along 
the [100] direction, the dashed lines and open circles in 



Fig. 4 demonstrate that for both point group symme- 
tries a [110] directed mode can occur. These sets of data 
support the idea of an interchange of the IGM vibra- 
tion direction as would be expected, for example, with 
hindered rotational motion of the excitation about the 
lattice site. This concomitant low frequency component 
of the IGM may provide new experimental ways to excite 
and identify these nonlinear excitations. 
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